#include <iostream>
#include <vector>
#include <string>
#include <time.h>
// #include <sstream>
#ifdef _MPI
#include <mpi.h>
#endif

#include "predefined.hpp"
#include "global.hpp"
#include "mesh.hpp"
#include "node.hpp"
#include "solution.hpp"
#include "flux.hpp"
#include "preprocess.hpp"
#include "solver.hpp"

using namespace std;

int main(int argc, char *argv[])
{
    #ifdef _MPI
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &iCPU);
    MPI_Comm_size(MPI_COMM_WORLD, &nCPU);
    #endif
    string caseFileName(argv[1]);

    int iStep, iStart, iRK, seRK, maxRK;
    double dt, tN, tRK, costTime, remainTime;
    time_t tStart, tEnd;
    string msg;
    Mesh msh;
    Solver solver;
    Preprocess preProcess;
    Node *node;
    ExNode *enode;
    BufferNode *bnode;
    Solution *sol;
    Flux *flux;
    #ifdef _LNS
    ViscousFlux *vflux;
    #endif
    Meanflow *mf;

    msg.reserve(MAXCHAR);
    logCtrl.setFileName(iCPU);
    logCtrl.welcome();
    msg = "This is the CPU core: " + to_string(iCPU) + " of " + to_string(nCPU);
    logCtrl.info(msg.c_str());
    msg = "The case control file is: "+caseFileName;
    logCtrl.info(msg.c_str());

    logCtrl.info("-------------- Reading Control File -------------");
    inputCtrl.readcaseCtrl(caseFileName);
    msg = "This is a " + to_string(inputCtrl.physDim) + " D problem.";
    logCtrl.info(msg.c_str());
    inputCtrl.setupParams();

    logCtrl.info("-------------- Reading Mesh -------------");
    msh.openFile();
    msh.readSize();
    msh.readBC();
    msh.subZone();
    msh.coreAssign(iCPU);
    msh.bounNodeCount();
    
    node  = new Node[msh.nLNode];
    sol   = new Solution[msh.nLNode];
    flux  = new Flux[msh.nLNode];
    #ifdef _LNS
    vflux = new ViscousFlux[msh.nLNode];
    #endif
    mf    = new Meanflow[msh.nLNode];
    msh.readCoor(iCPU, node);
    msh.nodeTyp(node);
    msh.closeFile();
    msh.readOversetMeshTab(msh);
    enode = new ExNode[msh.NExNode];
    msh.readOversetMeshInfo(iCPU, enode,msh);
    preProcess.calJacobi(msh, node);
    outputCtrl.writeJacobi(msh, node);
    preProcess.initialize(msh, node, sol, flux);
    solver.initializeBuffer(msh,node,bnode,sol);
    preProcess.initializeMeanflow(msh, node, mf);
    preProcess.meanflowCorrection(msh, node, mf);
    preProcess.spatialDiffMeanflow(msh, node, mf);
    outputCtrl.writeInitial(msh, node, sol,mf);

    dt = inputCtrl.dt;
    tN = 0;
    iStep = 0;
    iStart = 0;
    iRK = 0;
    if (inputCtrl.restartFlag==1) preProcess.readRestart(msh,sol);
    logCtrl.info("########################################");
    logCtrl.info("########## Entering Iteration ##########");
    for (iStep = iStart; iStep < inputCtrl.maxStep; iStep++)
    {
        if (iStep / inputCtrl.monitorFreq * inputCtrl.monitorFreq == iStep) time(&tStart);
        seRK = (iStep + 1) / 2 * 2 - (iStep); // 0  1
        #ifndef _RK46
        maxRK = 6 * seRK - 5 * (seRK - 1);    // 5  6
        #else
        maxRK = 6* seRK - 4 * (seRK - 1);    // 4  6
        #endif
        //maxRK=1;
        for (iRK = 0; iRK < maxRK; iRK++)
        {
            tRK = tN + solver.s.rk[seRK][2][iRK] * dt;
            solver.BCWallInviscidCorrectVelocity(msh, node,sol);
            solver.BCWallCorrectVelocity(msh, node,sol);
            solver.BCInflow(tRK, msh, node, sol, flux);
            solver.BCSymmetryPlane(msh, node,sol,flux);
            solver.BCNull(tRK, msh, node,sol,flux,mf);
            solver.spatialDiff(msh, node, sol, flux);
            solver.flux(msh, node, sol, flux, mf);
            solver.source(tRK, msh, node,flux);
            #ifdef _LNS
            solver.viscousFlux(msh, sol, flux, vflux,mf);
            #endif
            solver.BCWallInviscid(tRK, msh, node, sol, flux,mf);
            solver.BCWall(tRK, msh, node, sol, flux,mf);         
            solver.BCFarfield_Radiation(tRK, msh, node, sol, flux, mf);
            solver.artificialDamping(msh, node,sol,flux);
            solver.RK(iRK, seRK, msh, sol, flux);
            solver.exchSol(msh,node,enode,sol);
            solver.BCSymmetryPlane(msh, node,sol,flux);
        }
        if ((iStep - (inputCtrl.filterFreq - 1)) / (inputCtrl.filterFreq) * (inputCtrl.filterFreq) == (iStep - (inputCtrl.filterFreq - 1))){
             solver.filter(msh,node,enode,sol,flux);
            solver.BCSymmetryPlane(msh, node,sol,flux);
            solver.exchSol(msh,node,enode,sol);
            
        }

        if ((iStep - (inputCtrl.monitorFreq - 1)) / (inputCtrl.monitorFreq) * (inputCtrl.monitorFreq) == (iStep - (inputCtrl.monitorFreq - 1)))
        {
            solver.detectDivergence(msh, node, sol);
            time(&tEnd);
            costTime = difftime(tEnd, tStart) / 60.0;
            remainTime = (inputCtrl.maxStep - iStep) / inputCtrl.monitorFreq * costTime / 60.0;
            msg = "iStep: " + to_string(iStep+1) + " ct: " + to_string(tN) + " pt: " + to_string(tN * inputCtrl.time_ref) + " cost time (m): " + to_string(costTime) + " remain time (h): " + to_string(remainTime);
            logCtrl.info(msg.c_str());
        }
        solver.residual(msh,sol);
        tN += dt;
        if (outputCtrl.outTyp == 0)
        {
            if ((iStep>=outputCtrl.outSt) && ((iStep + 1) % outputCtrl.outFreq == 0))
            {
                msg = "Write solution " + to_string(outputCtrl.outIdx) + " time: " + to_string(tN);
                logCtrl.info(msg.c_str());
                outputCtrl.writeSol(iStep, tN, msh, node, sol);
            }
        }
    }
    logCtrl.info("Calculating finished.");
    delete[] node;
    delete[] enode;
    delete[] flux;
    delete[] sol;
#ifdef _MPI
    MPI_Finalize();
#endif
}